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ABSTRACT 

The cosmic reionization of hydrogen was the last major phase transition in the evolution of the 
universe, which drastically changed the ionization and thermal conditions in the cosmic gas. To the best 
of our knowledge today, this process was driven by the ultra-violet radiation from young, star-forming 
galaxies and from first quasars. We review the current observational constraints on cosmic reionization, 
as well as the dominant physical effects that control the ionization of intergalactic gas. We then focus on 
numerical modeling of this process with computer simulations. Over the past decade, significant progress 
has been made in solving the radiative transfer of ionizing photons from many sources through the highly 
inhomogeneous distribution of cosmic gas in the expanding universe. With modern simulations, we have 
finally converged on a general picture for the reionization process, but many unsolved problems still 
remain in this young and exciting field of numerical cosmology. 



Subject headings: cosmology: theory - large-scale structure of universe 
medium - methods: numerical - radiative transfer 



galaxies: formation - intergalactic 



1. Introduction 

Just as the 18th century explorers finished charting 
out most of the Globe, the 21st century explorers of the 
universe will most likely finish charting out the main 
areas of the cosmic evolutionary map. Already, the pre- 
vious century has witnessed some major discoveries, and 
over the past two decades we have started converging on 
a "Standard Cosmological Model" . Two main ideas lay 
in the foundation of modern cosmology: the expansion 
of the universe and the formation of cosmic structure. 
The expansion of the universe was first discovered by 
Edwin Hubble in 1929 and it was recently found that 
the exp ansion is accelerating due to a mysteriou s dark 
energy (jRiess et al.lll998l: iPerlmutter et al.lll999l) . This 
accelerated expansion continues to dilute the average 
density of matter in the universe, giving space an in- 
creasingly empty appearance. 

In contrast, structure formation involves gravitational 
contraction to higher densities. Tiny fluctuations in the 
density of matter at early times grew through gravi- 
tational instability to give rise to much larger cosmic 
structures at later times. The matter distribution, in 
which 80 — 85% is dark matter and 15 — 20% is cosmic 
gas, evolved to form a skeleton of high-density regions, 
called the "large-scale structure" or "cosmic web" . Em- 
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bedded in the highest density regions are the galaxies 
and clusters of galaxies that colour our picture of the 
cosmos. The large-scale distribution of galaxies in our 
cosmic neighbourhood has been cataloged by a number 
of surveys, culminating in the massive effort of the Sloan 
Digital Sky Survey (SDS^]). While some properties of 
the galaxy distribution are well understood within the 
framework of the Standard Cosmological Model, many 
questions remain about how galaxies form and evolve. 

As modern cosmologists strive to document the im- 
portant events in the cosmic history, they gain better 
understanding of how the two main ideas fundamen- 
tally shape the evolution of the universe. Currently, we 
have scoped out several periods of the cosmic evolution- 
ary map, but some major epochs have not yet been well 
charted. In this review, we focus on the exciting explo- 
ration of the epoch of reionization, a frontier in modern 
cosmology. 

In the cosmic chronology, the earliest observed infor- 
mation comes from the time when the universe was 
about 380,000 years old. At this stage, the cosmic 
plasma and radiation have cooled enough to allow elec- 
trons to combine with protons to form stable neutral 
hydrogen atoms. As a result, the primordial radiation 
mostly stopped scattering and since then has travelled 
largely unimpeded through the universe, to be detected 
by us as the cosmic microwave background. Measure- 
ments of the fluctuations or anisotropies in this radia- 
tion, first by the Cosmic Background Explorer (cobeH) 
satellite, followed by many other experiments, and re- 
cently to unprecedented precision with the Wilkinson 



■"-http:/ /www. sdss.org/ 
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Microwave Anisotropy Probe (WMAF(f|) satellite, have 
revealed much about the initial conditions of the early 
universe. 

After the process of recombination, the next several 
tens of millions of years are referred to as the "Dark 
Ages" . During this period, the matter in the universe 
continued to evolve under the influence of gravity, even- 
tually forming environments where early star-forming 
galaxies and rare quasar^ were born. These luminous 
objects produced ultra-violet photons that are energetic 
enough to dissociate the electron from neutral hydrogen 
atoms. As the ionizing radiation escaped into the space 
between galaxies, called the "intergalactic medium" , the 
reionization of the cosmic gas begins. 

The emergence of luminous sources marks the begin- 
ning of the epoch of reionization. A few hundred million 
years later, reionization is believed to be completed, for 
reasons which we discuss in more detail later. The cos- 
mic reionization of hydrogen was the last major phase 
transition in the evolution of the universe. During that 
epoch, helium was also ionized, but typically only one 
of two electrons are dissociated from each atom. The 
full ionization of helium is believed to have occurred 
at a later time when quasars, with photons energetic 
enough to dissociate the second electron, become suffi- 
ciently abundant. 

The study of reionization has emerged as a frontier 
topic in cosmology and is the focus of this review. We 
first discuss the current observational constraints in Sj2] 
and then describe the dominant physical effects that 
control the ionization of the intergalactic gas in fJ31 In 
2]we present a qualitative description of the numerical 
methods for modeling reionization and in |j5]wG summa- 
rize results from modern computer simulations. 

2. Observational Constraints on Cosmic Reion- 
ization 

Cosmic reionization by early galaxies and quasars left 
some residual neutral hydrogen and an abundance of 
free electrons in the intergalactic medium. These trac- 
ers have been used to study the epoch of reionization. 
The residual neutral hydrogen is probed through Lyman 
alpha absorption in the spectra of high redshift quasars, 
while the free electrons are detected through Thomson 
scattering of the cosmic microwave background. We now 
review the current major observational constraints and 
later in i j5.2.2l discuss how the neutral hydrogen can be 
observed directly through high-sensitivity radio obser- 
vations. 



^ http: / /lambda. gsfc. nasa. gov / product / wmap / 

^Quasars are active nuclei of galaxies, powered by supermassive 
black holes. While they are located inside galaxies, they are often 
treated as a separate class of sources, as the spectrum of ionizing 
radiation that they produce differs substantially from the typical 
stellar spectra of "normal" galaxies. 




Fig. 1 . — Measurements of the fraction of light from high 
redshift quasars, transmitted through the intergalactic 
medium; the rest of light is absorbed in the Lyman alpha 
line of neutral hydrogen. The rapid evolution at z > 5.5 
is commonly interpreted as the rapid evolution in the 
abundance of neutral hydrogen in the universe, possibly 
i ndicating the en d of the reionization epoch at z 6 
([Fan et al.ll2006a[) . Figure is courtesy of X. Fan. 



2.1. Lyman alpha absorption 

The first major constraint on cosmic reionization 
comes from observations of high-redshift quasars. The 
residual neutral hydrogen in the intergalactic medium 
can be quantified using Lyman alpha absorptioiH fea- 
tures in the observed spectra. Most of the very distant 
quasars have been discovered in the SDSS, although a 
few have been found by other searches. Wc refer the 



reader to a comprehensive review bv lFan et all (|2006al ) 
for a recent discussion of observational searches for high- 
redshift quasars. 

The absorption meas urements are summ arized in Fig. 
[TJ which we adopt from I Fan et al. ( 2006a ). At interme- 
diate redshifts (2 < z < 5) the residual neutral hydro- 
gen in the intergalactic medium causes absorption at the 
'--^ 10 — 50% level. The transmitted fraction gradually 
decreases at higher redshifts, but is still consistent with 
a highly ionized universe. However, the absorption in- 
creases much more rapidly between z = 5.5 and z — 6, 
and the few observations at z > 6 suggest perhaps even 
more rapid change. 

One commonly adopted interpretation of these obser- 
vations is that the universe was much more neutral at 
the epochs probed by the higher redshift data. The 



^ Since Lyman alpha absorption is resonant, the absorbed light 
eventually gets re-emitted; however, since it is re-emitted in a 
random direction away from the observing telescope, the scatter- 
ing process appears as absorption to an observer. We, therefore, 
adopt the widely used (although inexact) term "Lyman alpha ab- 
sorption" to describe this process. 
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transmitted fraction is observed to change by 2 — 3 
orders of magnitude in only a small redshift interval 
Az « 0.5. The rapidity of the change, if extrapolated 
further, would suggest that the reionization epoch ended 
not long before 2 = 6, perhaps somewhere between 
z — 6 and z = 6.3. This argument was presented shortly 
after t he first high-r edshift quasars were discovered by 



tional support fe.g. iFan et al 



SDSS ( Becker et al.l 2001i). and has since drawn addi- 



2OO2I: I White et al.ll2003l: 



Gnedinll2004l: iFan et al.ll2006bHGnedin k Fanll2006l) . 

However, the epoch of reionization is not directly 
probed by these observations. Since the Lyman al- 
pha opacity of the fully neutral intergalactic medium 
at these redshifts is very large, of the order of 10^ or 
even higher, it only takes a small neutral fraction to 
have nearly complete absorption. The lowest measured 
transmitted fraction of ~ 10""^ at 2: « 6 is consistent 
with a volur ne-averaged neutral fraction of only 10"'* 
to 10-3 (e.g. iLi dz et al]l2006l : [Becker et aLlbOOTT ). The 
evolution of the neutral fraction at higher redshifts and 
the exact timing of reionization remains, at present, un- 
known. 

2.2. Thomson electron scattering 

The second major constraint comes from observations 
of the cosmic microwave background radiation. As the 
universe evolves from the recombination epoch to the 
present time, the cosmic microwave background pho- 
tons Thomson scatter with electrons dissociated dur- 
ing the reionization epoch. The scattering results in 
a small suppression of cosmic microwave background 
anisotropics on all scales and also generates additional 
polarization on large angular scales. A convenient quan- 
tity that parametrizes both effects is the total optical 
depth, tt, to Thomson scattering. 

The WMAP satellite has detected the effects of Thom- 
son scattering due to reionization and the first measure- 
ment of Tt came after only one year of observations. Th e 
value reported, tt = 0.17 ± 0.04 (jKogut et al.ll2003h . 
was unexpectedly large. Considering that the contri- 
bution to Tt of the fully ionized intergalactic medium 
between z = and z = 6 is only 0.04, a surprisingly 
large contribution, 0.13 ± 0.04, was left for the reion- 
ization epoch. Such a large value for tt would indicate 
an unusually prolonged or very early reionization epoch 
that is difficult to understand within the framework of 
the Standard Cosmological Model. 

However, recent measurements by WMAP after more 
years of observations are significantly lower with smaller 
uncertainties. The latest value, tt = 0.087 ± 0.017, is 
based on 5 years of data, and is in good agreement with 
the models of reionization that are based on the modern 
state-of-the-art numerical simulations, as we discuss in 

m 

With a few more years of observations, WMAP will 
have slightly improved measurements of tt. Further- 
more, a much more precise measurement by the upcom- 



ing Planck Surveyor missior|§ (to be launched in 2009) 
will allow the Thomson optical depth to be used to ro- 
bustly constrain models of reionization. 

2.3. Lyman alpha emitters 

The third observational probe of reionization is pro- 
vided by a class of high-redshift galaxies known as Ly- 
man alpha emitters. Star-forming galaxies emit a signif- 
icant fraction of their radiation at Lyman alpha wave- 
lengths and these galaxies are recognized by the strong 
Lyman alpha emission lines. During the reionization 
epoch, the Lyman alpha emission from these galaxies is 
scattered by surrounding neutral hydrogen. The scat- 
tering changes the distribution of apparent brightness 
for the population of Lyman alpha emitters. 

Changes in the abundance and distribution of Ly- 
man alpha emitters can be used to detect the tran- 
sition from the largely ne utral to largely ionized in- 



tergalactic medium (e.g Malh otra fc Rhoadsl 20ol: 



Kashikawa et al.ll2006l : IStark et al...2007 : .McQuinn et al 



2007al) . However, recent research has shown that the 
interpretation of those measurements is extremely com- 
plex, and no consensus exists at present on how the 
observations of Lyman alpha emitters should be used 
to constrain reionization. Additional observations and 
better theoretical understanding of t he effects of scat- 
tering on Lyman al pha emitters (e.g. Tasitsiomi 2006t 
Diikstra et al.ll2007( ) are required to place more reliable 
constraints on the reionization epoch. 

3. Physics of Reionization 

In modeling reionization, the main important physical 
process is the transfer of ionizing radiation through the 
inhomogeneous distribution of cosmic gas. Since ion- 
izing photons are not scattered, but only emitted and 
absorbed, and travel in straight lines between the acts 
of emission and absorption (with a minor and often neg- 
ligible complication that they redshift as they travel in 
the expanding universe), the physics of reionization pro- 
cess is largely determined by the physical properties of 
sources and sinks for the ionizing radiation. 

3.1. Sources of ionizing radiation 

The nature of reionization sources remains poorly con- 
strained. Two main types of astrophysical sources, mas- 
sive stars in star- forming galaxies and quasars, are of- 
ten considered likely to be the dominant ones. How- 
ever, other more exotic possibilities, like decaying dark 
matter particles or evaporating primordial black holes, 
cannot be completely discounted yet. We will restrict 
our focus hereafter to galaxies and quasars as the only 
two choices of reionization sources since all simulations 
of reionization to date have only considered these cases. 
The reader must be cautioned, however, that the reality 



' ]http: //www.rssd. esa. Int/lndex .php?project=planck| 
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Fig. 2. — Plausible ranges for the number of ioniza- 
tions per hydrogen atom from quasars (red) and galax- 
ies (blue) - the detailed explanations of adopted limits 
are given in the text. The large black dot is a necessary 
reionization conditions of one ionization per hydrogen 
atom by z = 6. 



can be more complex than reionization modelers have 
been willing to accept so far. 

The most basic condition on any type of source is 
that the total number of ionizations produced by these 
sources must at least equal the number of hydrogen 
atoms in the universe. Otherwise, it will not be possible 
for that kind of source to reionize the whole universe. 
This condition is actually only a minimum requirement. 
Since the cosmic gas is able to recombine, more than 
one ionizing photon per hydrogen atom is necessary to 
maintain the ionization. The required photon to atom 
ratio remains a major unknown - computing this quan- 
tity is an important task for the computer simulations 
that we describe below. 

Fig. [2] shows a plausible constraint on the number of 
ionizations by galaxies and quasars. Since we lack a 
comprehensive theory of reionization, we show for each 
type of source a range of estimates that most likely 
bound the truth. The lowest quasar estimate is based 
on the direct extr apolation of the obse rved quasar lu- 
minosity function ([Hopkins et al.l 120071 ) to higher red- 
shifts, assuming that each ionizing photon produces a 
single ionization. Quasars are believed to emit all of the 
ionizing radiation they produce and therefore, the total 
ionizing luminosity can be easily estimated by knowing 
their abundance. However, this extrapolation is likely 
to be an underestimate for two reasons. First, the spec- 
trum of ionizing radiation from quasars is quite hard, 
with a large number of energetic photons and a mean 
photon energy of about one keV. When such an en- 
ergetic photon ionizes a hydrogen atom, an energetic 
electron is produced. These energetic electrons can, in 
turn, ionize more atoms. Thus, a single ionizing pho- 
ton can ionize more than one atom via the "secondary 



ionizations" process. While the exact number of sec- 
ondary ionizations requires complex modeling, a rule of 
thumb for a crude estimate is that one ionization is pro- 
duce d for every 40 eV of th e ionizing photon energy 
(,Shull fc van Steenber j|l985l ). 

Second, observations miss the low-luminosity objects 
that are fainter than the detection limits. In addition, 
the rarest bright ones may not be found within the fi- 
nite surveyed volumes. Thus, the observations are likely 
to underestimate the total number of ionizing photons 
coming from all quasars. While the exact magnitude of 
this underestimate is difficult to compute, simple mod- 
els indicate that it is unlikely to be more than a factor of 
3 to 5 (M. Volontcri, private communication). The up- 
per bound for the quasar range is therefore obtained by 
simply multiplying the lower bound by a factor of 100, 
and it should be considered an extremely conservative 
upper limit. Thus, Fig. [2] illustrates a co nclusion that 



is well known since the pioneering work of iMadau et al 



(1999); quasars alone are not powerful enough to reion- 
ize the whole universe by z = 6, as required by the 
observations discussed in ^ 

With star-forming galaxies the situation is more com- 
plex. On one hand, the luminosity functions of galax- 
ies are measured with reasonable precision all the way 
to z = 6, and some measurements exist up to z ~ 9 
( Bouwens et al. 20081) . The uncertainty due to incom- 



plete observations is, therefore, not expected to be large 
(Gnedin 2008 ). On the other hand, galaxies do not emit 
all of the ionizing radiation they produce and the total 
ionizing luminosity can not be estimated by just know- 
ing their abundance, unlike for quasars. Galaxies only 
leak a fraction of the radiation produced by their stars, 
with the rest being absorbed locally by gas clouds in the 
interstellar medium. 

The fraction of ionization radiation that escapes from 
galaxies, /esc, is a major uncertainty. The upper bound 
in Fig.[2]is computed from the observed galaxy luminos- 
ity functions with the assumption that /esc — 20%. This 
assumption may not be reasonable, however, since ob- 
servations of bright galaxies at lower redshifts indicate 
that /esc ~ 2%. Thus, the upper bound for the stellar 
contribution to the ionizing budget should be, just as 
one for quasars, treated as an extremely conservative 
upper limit. The bottom bound for the stellar contri- 



bution follows the model of iGnedin et al.l (j2008l ). which 
is based on the combination of observational data and 
numerical simulatio ns of high red s hift galaxies. If the 



simulations used bv I Gnedin et al. ( 2008f ) are not accu- 



rate enough, the lower bound may be an underestimate, 
although it is unlikely to be off by more than a factor 
of 3. 

A simple conclusion can be drawn from Fig. [2l Stars 
in galaxies are an important and most likely domi- 
nant source of ionizing radiation during the reionization 
epoch. However, if the escape fractions from galaxies 
are as low as is indicated by numerical simulations, then 
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Fig. 3. — Correlation functions of galaxies (left) and 
dark matter (right) as a function of the comoving dis- 
tance r at four different epochs from the Millennium 
simulation of structure formation (jSpringel et al. I l2006h . 
Figure is reprinted with permission by V. Springel. 

quasars can contribute as much as ^ 50% of all ioniza- 
tions of hydrogen. 

To complicate things even further, galaxies and 
quasars are not distributed randomly in the universe, 
but populate a complex large-scale structure of clusters, 
groups, filaments, and voidq3- The galaxies and quasars 
are clustered and aggregate together, such that the typi- 
cal separation between nearest neighbors is smaller than 
what is expected if they are randomly distribute41- The 
simplest measure of any clustered distribution is a "cor- 
relation function", ^(r), which describes the fractional 
excess of neighbors at a given distance r over a random 
distribution. In Fig. [3] we show the galaxy correlation 
function at a range of redshifts as well as the correla- 
tion function for the dark matter at the same redshifts. 
These two functions are not the same as galaxies are al- 
ways more clustered than the dark matter. The galaxies 
are said to be "biased" and at high redshifts, the bias- 
ing can be quite large, as Fig. [3] illustrates. In §3.31 we 
discuss the effects of the clustering of sources on the 
reionization process. 

3.2. Sinks of ionizing radiation 

During the epoch of reionization, each neutral hydro- 
gen atom serves as a sink for an ionizing photon, or at 
least a fraction of an ionizing photon if secondary ion- 
izations are taken into account. If this was the only 
type of sink, the theory of reionization would have been 
completed by now. 

Unfortunately for theorists, there exist other types of 
sink for ionizing photons. For example, the already ion- 



Qualitatively and even quantitatively, the distribution of matter in 
the universe is similar to the distribution of human population on 
the Earth. There are huge metropolitan areas with the population 
density orders of magnitude higher than the mean, large and small 
cities and towns strewn mostly along major highways, and small 
villages and hamlets in the most sparsely populated areas. 
*If all the people on Earth were randomly distributed, the average 
distance to your nearest neighbor would be 160 m, but most likely 
you will find someone much closer if you look around right now. 



ized gas can recombine if its density is high enough. 
These recombinations will serve as an additional sink 
of ionizing photons, since each newly recombined atom 
would need to be ionized again to complete the reioniza- 
tion of the universe. It is generally believed that recom- 
binations can consume a measurable but not the dom- 
inant fraction of all ionizing photons. They are an im- 
portant sink for ionizing radiation at the earliest stages 
of reionization, but become progressively less important 
as the universe expands and the gas density decreases. 

Another type of sink for ionizing radiation, called "Ly- 
man limit systems" , exists both during and after the 
reionization epoch. Lyman limit systems are observed 
in the spectra of distant quasars as absorption systems 
with column densities of neutral hydrogen exceeding 
about lO^^cm"^. At these column densities, the optical 
depth for ionizing radiation is t^l ^ 1 at the Lyman 
limit wavelength A = 912A. 

Usually, it is extremely difficult to make a connection 
between absorbing systems in quasar spectra and the 
physical objects in space that correspond to them. At 
present, astronomers know some properties of Lyman 
limit systems, for exam ple, that they are mostly ionized 
and not highly neutral ( Prochaska 19991 ). However, we 
still lack the knowledge about what kind of physical ob- 
jects they are, or how they correlate with galaxies and 
quasars. 

For studies of reionization the most important factor 
is the abundance of Lyman limit system as a function 
of redshift. After all, the only property of sinks of ion- 
izing radiation is that they are indeed sinks, and that 
they absorb all ionizing photons that hit them, and it 
is not that important whether they are clouds of hydro- 
gen gas or photon-eating space-dwelling monsters. For 
example, it generally makes little difference whether the 
Lyman limit systems are highly clustered or distributed 
unformly in space, as long as they have the same abun- 
dance per unit redshift. The abundance mainly deter- 
mines the probability for an ionizing photon to be ab- 
sorbed, irrespectively of the true nature or spatial dis- 
tribution of the absorbers. 

The measurements of the abundance of the Ly- 
man limit systems at z < 4 has been presented by 



Storrie-Lombardi et al.l (| 19941 ) some time ago, although 



more recent and accurate measurements from the SDSS 
should be forthcoming shortly. However, a more physi- 
cally intuitive quantity is the "mean free path" of ion- 
izing radiation Amfp, defined as the average distance 
an ionizing photon travels through space before being 
absorbed. The measured abundance of Lyman limit 
systems can be easily converted into the measurement 
of the mean free path of ionizing radiation. At z ~ 4, 
Amfp « 90h-^ Mpc in comovingj units, with about 15% 



''Comoving distances are defined as real, physical distances scaled 
by 1 + 2 to account for the expansion of the universe; they are con- 
venient to use because comoving distances remain constant with 
time in an inertial reference frame. In addition, in order to scale 
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uncertainty (l Miralda-Escudell2003h . 

At higher redshifts the constraints do not yet exist, 
but the observed evolution of the Lyman limit abun- 
dance is well described by a power- law in (1 -f z), so 
it can be directly extrapolated to higher redshifts for a 
plausible, but uncertain, estimate of the mean free path 
from the Lyman limit systems at earlier times. 

3.3. General overview of reionization process 

Armed with basic knowledge of sources and sinks for 
ionizing radiation, we can try to understand general fea- 
tures of the reionization process. While our goal is not 
to give a comprehensive review of the theory of reion- 
ization, it is useful to understand the main process and 
relevant spatial scales, as these are important for un- 
derstanding the limitations and the degree of scientific 
fidelity of simulations that we describe below. 

The story of reionization is not difficult to imagine. 
As sources of radiation begin to produce energetic pho- 
tons, they will start ionizing regions around them, usu- 
ally called "HII regions" (after the spectroscopic nota- 
tion for ionizing hydrogen, HII) or, more colloquially, 
"ionized bubbles" . The ionized bubbles keep expanding 
until... well, until that expansion stops. 

There could be two physical reasons why the expan- 
sion stops: either because all the ionized bubbles over- 
lap, and all of the intergalactic medium gets reion- 
ized, or because the Lyman limit systems inside the 
bubbles absorb all ionizing photons available for ion- 
izing gas outside a bubble. Since the bubbles originate 
around sources, the choice between the two scenarios is 
determined by the relationship between the mean free 
path for ionizing radiation and the average distance be- 
tween ionizing sources. If the source separation is much 
smaller than the mean free path ("abundant sources" 
scenario), then most of ionizing bubbles overlap before 
the absorptions by Limit limit systems becomes signifi- 
cant, the overlap of bubbles happens fas t, and results in 
a quick reionization of the universe g.e. lGnedinlbnnnh . 

If the opposite is true, and the sources are so 
strongly clustered that their correlation length is much 
larger than the mean free path ("rare sources" sce- 
nario), then ionized bubbles reach their maximum sizes 
(comparable to a few times the mean free path) be- 
fore they can overlap , and r ei onization process stalls 



( Miralda-Escude et al. I2OOOI : Furlanetto fc Mesingeil 



20091 ). The whole universe can then be reionized only 
when either the abundance of Lyman limit systems 
decreases (due to the cosmic expansion and the associ- 
ated density decrease or because they get ionized more) 
or sources become more numerous or less strongly clus- 
tered. In other words, reionization does not actually get 
completed in a pure "rare source" scenario. Only when 



out the dependence of cosmic distance on the insufficiently accu- 
rately known Hubble constant Hq, the comoving distances have 
been customarily expressed for many years in Mpc, where 
h = _H"o/(100km/s/Mpc). 




Fig. 4. — The estimate of the mean free path from 
the Lyman limit systems as an extrapolation of the 
observed evolution of the Lyman limit abundance 
( Storrie-Lombardi et al. 1994 , ; red band). The blue 



band shows the mean galaxy sep a ration as a function 



of redshift (|Bouwens et all [20071 . l2008f) . Black dots 



with error-bars are observational estimates of the mean 



free pa th from the spectra of SDSS quasars ()Fan et al 



2006b|). The dashed red line shows an extrapolation of 



the mean free path from the Lyman limit system that 
matches the SDSS estimates. 



sources become sufficiently abundant, with at least one 
source (or a group of clustered sources) per every sphere 
with radius of a few Amfp, can reionization of the whole 
universe end. 

In principle, these two scenarios can be distinguished 
by comparing the abundance of Lyman limit system (i.e. 
the mean free path) and the abundance of sources dur- 
ing reionization. In practice, the situation is rather com- 
plex. The current observational constraints are summa- 
rized in Fig. 2) There we show the observed mean free 
path from Lyman limit systems (red square at z = 4) 
and its extrapolation to higher redshifts (red hatched 
band) using the observed rate of evolution of Lyman 



limit systems at lower redshifts (jStorrie-Lombardi et al 



11994) . the estimate of the mean separation between 
galaxies from recent survey s of Lyman Break galaxies 
( Bouwens et al.ll2007ll2008f) . as well as observational es- 
timates of the mean free path fo r ionizing radiatio n from 
the spectra of SDSS 
points with error-bars 

Taken at face value, the existing constraints seem to 
suggest that the mean free path is larger than the mean 
galaxy separation, and that the "abundant sources" sce- 
nario is realized. For example, the fact that the estimate 



uasars ( Fan et al 
lol 



2006bl . black 



^"The red hatched area is computed assuming dN/dz = (3.3 ± 
0.5) ((1 -I- z)/b)'^-^^^°-^^; the blue hatched area is obtained from 
iBouwens et al. fits as 4>* the dashed red line as- 

sumed the Lyman limit system evolution in the form dN/dz = 
3.3{(l + 2)/5)^ 
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of the mean free path from the spectra of SDSS quasars 
drops rapidly at z « 6 has been used to argue that the 
fast overlap of ionized bubbles occurred at that redshift 
(iGnedin k Fanll2006l) . 

The situation, however, is likely to be more complex. 
First of all, both the mean free path Amfp and the mean 
galaxy separation dsRC should be taken as characteristic 
scales, not as exact distances. Since the mean galaxy 
separation is only a factor of few below the extrapolated 
mean free path, it is not really in the regime where it is 
"much smaller" than Amfp- 

In addition, galaxies are highly clustered, as Fig. [3] 
illustrates. In a highly clustered distribution, the mean 
distance between galaxies is a poor indicator of actual 
spacing between sources: some galaxies will be so close 
to each other that they effectively form a single, more 
powerful source, while clusters of sources can be spaced 
by much more than dsRC- For example, in the present 
day universe the average distance between galaxies is 
about 4 — 5h~^ Mpc; never-the-less, voids in the large- 
scale distribution of galaxies reach sizes of several tens of 
Mpc and filamentary super-clusters approach hundreds 
of Mpc in length. 

It is also possible that our extrapolation of the mean 
free path from the Lyman limit systems is incorrect. For 
example, if the abundance of the Lyman limit systems is 
larger at higher redshifts, the mean free path would be 
smaller. As an illustration, the dashed red line in Fig. 
|4] shows the extrapolation of the Lyman limit system 
abundance that matches well the estimate of the mean 
free path from the SDSS quasars. In that case the mean 
free path gets smaller than the galaxy correlation length 
at z « 6.5. 

Thus, we must conclude that the existing knowledge 
of the reionization epoch does not clearly prefer either 
abundant or rare sources scenario; the reality is, most 
likely, to be somewhere in between the two extremes. 
That makes the theorist's task of creating simple, an- 
alytical theories of reionization so much harder, and 
points toward detailed numerical simulations of reion- 
ization as the ultimate method of choice. 

4. Numerical Methods 

Computer simulations provide a powerful and versa- 
tile tool for solving the fundamental physics of gravi- 
tation, gas dynamics, and radiative transfer, which are 
important for modelling reionization. Many specialized 
algorithms have been developed and applied to simulat- 
ing the complex interactions between dark matter, gas, 
stars, and radiation. Numerical modeling has also ben- 
efitted from the rapid development in supercomputing 
technologies. Computing capabilities, such as proces- 
sor speed and memory capacity, have been growing ac- 
cording to Moore's Law, doubling approximately every 
two years. Together, these advancements enable increas- 
ingly larger and more realistic simulations to be run. 

With modern cosmological simulations, we can now 



study how the small, initial perturbations in the 
dark matter and cosmic gas distributions undergone 
nonlinear gravitational collapse t o form the large- 
scale structure of t he universe (see Bertschingei 1998t 
Springel et al. 2006i for reviews). The large-scale struc- 
ture, characterized by filaments and halos, sets up cos- 
mic environments where the early stars and galaxies 
form. While astrophysical processes such as star for- 
mation still can not be simulated from first principles, 
simulations do allow a more straightforward implemen- 
tation of intricate prescriptions, motivated by theory 
and calibrated against available observations. 

Over the past two decades, there have been consid- 
erable interest in cosmological simulations with radia- 
tive transfer. Earlier work focussed on solving for the 
radiation field around a single point source, which is 
applicable to the first stars or first galaxies. In the 
last few years, the attention has shifted towards under- 
standing how the larger distribution of early stars and 
galaxies photo-ionized and photo-heated the intergalac- 
tic medium during the epoch of reionization. Radiative 
transfer in cosmological simulations is still in its infancy 
compared to dark matter and gas dynamics, but it is 
rapidly gaining strength and scope. In the following 
sections, we review the computational methodology and 
requirements for simulating cosmic reionization. 

4.1. N-body and Hydrodynamic Algorithms 

Currently, there are two classes of cosmological sim- 
ulations. One class, designated "N-body" simulations, 
models the evolution of all matter in the universe as 
a collisionless fluid influenced only by gravitational dy- 
namics. The other class, called "N-body -I- hydro" or 
just "hydro" simulations, models both the collisionless 
dynamics of dark matter and the coUisional dynamics of 
cosmic gas. 

N-body algorithms are normally used to evolve the 
dark matter distribution, which is discretized into N 
number of particles of fixed mass, each with known po- 
sition and velocity. Newton's equations of motion can be 
solved in discrete time steps when given the acceleration. 
The exact force on a given particle can be calculated by 
doing a direct summation of the pair- wise forces exerted 
by all other particles. However, this is prohibitively ex- 
pensive for cosmological simulations because the num- 
ber of calculations scales as 0{N'^) wh ere N can be 



large , currently up to 7 billion particles (jTevssier et al 



20091 : Kim et al. 2008 ). Fortunately, several successful 
techniques have now been implemented to solve Pois- 
son's equat ion for gravity wit h reasonable 0{N log N) 
scaling (see iBertschinger 19981 for a review). 

Hydrodynamic algorithms solve the fluid equations 
for the cosmic gas, using grid-based ("Eulerian") or 
particle-based ( "Lagrangian" ) techniques. In the Eule- 
rian approach, the conservation equations for gas mass, 
momentum, and total energy are solved on a structured 
or unstructured grid of cells. The mass in a cell can have 
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practically any value and hence grid-based algorithms 
are said to have high dynamic range in mass. Eulerian 
codes using uniform grids are usually more suited to sim- 
ulating the intergalactic medium, which contains both 
dense and underdense gas. However, uniform grids have 
limited spatial dynamic range. Adaptive mesh refine- 
ment codes additionally have the capability to resolve 
small-scale structure like gas in halos. 

In the Lagrangian approach, fluid elements are rep- 
resented by particles of fixed mass and smooth particle 
hydro dynamics fSPH dLucvlll977 : Gingold fc Monaghan 
solvers are used to follow their trajectories. For 
each collisional particle, the dynamical forms of the 
fluid equations for density, velocity, and temperature are 
solved while smoothing over a small number of neigh- 
boring particles. SPH is especially suited for simulating 
small-scale structure because the Lagrangian flow natu- 
rally allows for high spatial dynamic range in high den- 
sity regions. On the other hand, underdense regions in 
the intergalactic medium are resolved with fewer parti- 
cles. 

Hydrodynamic algorithms can be run simultaneously 
with N-body algorithms to couple the gas and dark 
matter through their mutual gravitational attraction. 
Cosmological simulations of reionization also solve for 
atomic processes such as ionization, recombination, 
cooling, and heating. The photo-ionization and photo- 
heating of the gas are straightforward to calculate given 
the inhomogeneous radiation field. 

4.2. Radiative Transfer Algorithms 

Radiative transfer algorithms solve the evolution of 
the radiation field, taking into account emission, absorp- 
tion, and scattering processes. In general, the evolution 
is described by a differential equation for the specific 
intensity, which is a function of seven variables: 3D po- 
sition, 2D angular coordinates, time, and frequency. Be- 
cause of the high dimensionality of the problem, direct 
numerical solutions are computationally difficult and ex- 
pensive. 

For a discrete number TYrt of radiative transfer reso- 
lution elements, a direct solution requires 0{N^!^) oper- 
ations per frequency bin per time step. With this costly 
scaling, only low resolution simulations can be run if a 
brute-force solution is attempted. In order to be feasible 
for use in cosmological simulations, radiative transfer al- 
gorithms should scale close to linearly with the number 
of resolution elements, just like good N-body and hydro- 
dynamic algorithms. However, satisfying this criterion 
requires some level of physical approximations and com- 
putational optimizations. 

Existing algorithms can be broadly divided into three 
categories: moments , Monte C a rlo, an d ray-tracing 
methods. Recently, llliev et al] ( 2006a ) conducted a 



some results from a classical test of an ionized bubble 
expanding from a single source into a gas of uniform 
density and temperature. Below, we provide a qualita- 
tive description of the three main methods for modeling 
cosmological radiative transfer. 

4-2.1. Moments methods 

The radiative transfer equation for the specific inten- 
sity can be simplified by considering moments of the 
radiation field. It can be reduced to a simpler system 
of conservation equations for the photon energy den- 
sity and fiux. This is analagous to replacing the Boltz- 
mann equation for the fiuid distribution function with 
the Euler conservation equations for gas mass, momen- 
tum, and total energy. Furthermore, similar to the Euler 
equations having source terms on the right-hand side, 
for example coming from gravity, the radiative trans- 
fer moments equations have three important terms, two 
of which come from the emission and absorption. The 
third term, called the Eddington tensor, is related to the 
radiation pressure and is necessary to close the system 
of partial differential equations. 

Radiative transfer moments algorithms are naturally 
coupled to Eulerian hydrodynamic codes. Since the ra- 
diative transfer moments equations resemble the Euler 
hydro equations in form, they can be solved using the 
same, well-established techniques. Algorithms generally 
scale as C'(A'fit fog A^rt), independently of the number 
of sources. The source function can easily incorporate 
both point-source and diffuse radiation. This advanta- 
g eous feature is n o t gen erally shared by other methods. 

Gnedin fc Abell (1200 ih were the first to develop an al- 



gorithm specifically for reionization and proposed the 
Optfoally Thin Variabfo Eddington Tensor (OTVET) 
approximation. Because computing the exact Edding- 
ton tensor is a highly nontrivial task, the OTVET ap- 
proximation computes the tensor under the assump- 
tion that the absorptions are negligible, but then used 
this "optically thin" tensor in the full, "optically thick" 
equation for the radiation energy density and flux. That 
ensures the conservation of photon number and fiux, 
but causes errors in the direction in which the radiation 
fiux is advected. One positive is that these errors re- 
main under control at all times, and the accuracy of the 
OTVET approximation can always be estimated; if this 
accuracy is found to be inadequate, a different, more 
precise method must be used instead. Fig. [5] demon- 
strates that, overall, the OTVET approximation agrees 
well with other techniques, but, just like Monte Carlo 
methods, is somewhat more diffusive than direct ray- 
tracing algorithms. 

Several alternative schemes for computing the Ed- 
dington tensor have been developed recently (e.g. 



comparison of 11 cosmological radiative transfer codes 
using 5 simple test problems and found good general 
agreement between different algorithms. Figure [5] shows 



Aubert fc Tevssiej2008l:lFinlator et al.l2009l:|Petkova fc Springe! 



20081 ) . These various algorithms mainly differ in how the 
Eddington tensor is computed, and what assumptions 
are made in that computation. 
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Fig. 5. — A comparison of different numerical approaches for modeling radiative transfer bv llliev et al.l ( 2006al ). Each 
panel shows an octant of a spherical ionization bubble blown up by a single source, with the approach/code used in 
the simulation labeled above. For this classical test, the ray-tracing algorithm (C^-RAY) provides the most accurate 
solution; the moment method (OTVET) is usually somewhat diffusive, while the Monte-Carlo approach (CRASH) 
introduces random fluctuations that break the perfect symmetry of the problem. Figure is courtesy of I. Iliev. 



4-2.2. Monte Carlo methods 

In Monte Carlo methods, the radiative transfer equa- 
tion is solved using a probabilistic technique where the 
values of variables are randomly sampled from known 
probability distribution functions. The radiation field is 
discretized using photon packets and a number of pack- 
ets is emitted from each source. For each packet, initial 
conditions such as the emission location, propagation 
direction, and frequency are determined by randomly 
sampling from the appropriate probability distribution. 
As a packet is transported away along a radial path, it 
will intersect gas elements and for each crossing, a frac- 
tion of the photons is consumed based on the probability 
for absorption. 

For a single source, the number A^p of packets emitted 
should be comparable to the number Nrt of radiative 
transfer resolution elements in order to have informa- 
tion about the radiation field everywhere. For a gen- 
eral problem with Ns sources, the total number of pack- 
ets emitted is then NpNs. In reionization simulations, 
Ns can be very large, scaling linearly with N^t, and 
thus making the overall scaling be 0{N^r^). In prac- 
tice, optimizations are introduced such that the total 
number of packets is proportional to the number of res- 
olution elements and only 0{Njn) operations are done. 
This is usually achieved by reducing Np and comes at 
the expense of degrading the angular resolution around 
sources. Convergence tests, where A'^p is varied for ex- 
ample, must then be conducted to ensure that statisti- 
cal fluctuations do not significantly affect results. With 
good optimization, Monte Carlo methods can be effi- 
cient and reliable. 



CRASH (jCiardi et al.l 1200X1 : iMaselli et al.l l2003l ) was 
the first Monte Carlo code written for reionization. 
The latest version operates on static, grid-based den- 
sity fields and solves for the time evolution of H and He 
ionization fractions and gas temperatures. When com- 
pared with other cosmological radiative transfer codes. 



it showed good overall agreement on test problems 



but tended to have thicker ionization fronts (Iliev et al 



2006al ). In Fig.jSjthc broadening of the ionization front 
comes from the probabilistic nature where there is al- 
ways a finite chance that a photon packet stops before or 
travels beyond the actual ionization front. More recent 
imple mentations (e.g. Semclin ct al. 20G'3; lAltav et al 



20081) have been designed to couple to SPH, allowing 
better spatial resolution in high density regions. 

4.2.3. Ray-tracing methods 

Ray-tracing is the most popular of the three categories 
and there is a rich diversity of approaches. We first re- 
view the basic techniques and then discuss adaptive im- 
provements. For the basic case, we consider ray-tracing 
from a single source through a radiative transfer grid 
(e.g. lAbel et al.|[l999f) . A fixed number of rays is gener- 
ated with an isotropic distribution of propagation direc- 
tions. Walking downstream away from the source, each 
ray is cast into segments as it traverses the radiative 
transfer grid, basically one segment for every cell inter- 
sected. Photons are consumed in each segment based 
on the local optical depth for photo-ionization, which 
depends on the length of the segment and the density 
of absorbers in the cell. Each ray is calculated indepen- 
dently, and gets terminated when its photon count goes 
to zero. 

In a multi-source problem, the situation is more com- 
plicated because rays cannot be calculated indepen- 
dently in general. For a causal solution, all rays must 
be synchronized and traced in discrete steps, by one 
cell length at a time over a time step equal to the 
time it takes light to cross a cell. This gives the 
most accurate results, but is more costly as many 
steps are needed because the light-crossing time can be 
quite small compared to the duration of reionization. 
However, practical but approximate solutions do exist 
where sources can be processed independently, allowing 
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for much longer time steps (e.g. ISokasian et al.l 12001 
iMellema et~al]|2006al) . 

Ray-tracing is computationally expensive and also 
runs into the 0(7Vprp) scaling problem we discussed ear- 
lier for Monte Carlo methods. In the brute-force ap- 
proach, the number of rays cast per source must be 
equal to or greater than the number -/Vrt of radiative 
transfer resolution elements. This is necessary to have 
rays intersect cells far from the source, but it results 
in much more rays than necessary nearer to the source. 
Fortunately, there are adaptive improvements to make 
algorithms run faster and scale better. 



Abel fc Wandeltl (|2002l ) suggested an adaptive ray- 
splitting technique where the angular resolution contin- 
ually increases farther away from sources. For a given 
source, a small number of parent rays are cast and 
travel a short distance before they split into 4 daughter 
rays. Successive generations of splitting continue down- 
stream. Furthermore, adaptive ray-tracing algorithms 
can be made to scale more linearly with the number 
of resolution elements by gr ouping neighboring sources 
(iRazoumov fc Car dall' '2005*) . imiting the splitting of 



rays (IMcQuinn et al.. .2007b) . or merging near-parallel 



rays (|Trac fc Cenll2007l ). 



So far we have only discussed ray-tracing techniques 
for grid-based hydrodynamic codes. Several ray-tracing 
algorithms have also been developed for c osmological 
SPH simulati ons (e.g. lAlvarez et all l2006t ISusal l2006t 
Pawlik fc Schave 2008). The ray casting techniques in 



SPH are different because of the irregular distribution 
of the particles. One common approach is to emit a 
fix number of rays from each source and collect them 
using an unstructured grid based on the particle distri- 
bution. High spatial dynamic range is achieved for both 
the particle s and the rays in high density regions. 

C^-RAY (jMellema et al.l l2006al ) is one example of a 
grid-based ray-tracing code that has been used for simu- 
lating reionization. This photon-conserving code traces 
rays away from every source to every radiative trans- 
fer cell. While each source is processed independently 
and in random order, an iterative procedure is used to 
converge to a causal and accurate solution. The cur- 
rent version scales linearly with the number of sources, 
which is costly when Ng scales with A^rt- In Fig- [3 the 
ionized bubble is in good agreement with the analytical 
solution and with results from other ray-tracing codes 
(llliev et al. 2006a). 



4.3. Physical Scales and Computational Re- 
quirements 

Computer simulations of reionization require a large 
dynamic range in both length and mass scales. Even 
with modern supercomputers, only a portion of the en- 
tire range of scales can be probed by any single sim- 
ulation. High resolution is required to resolve small- 
scale structure such as radiation sources and sinks, the 
two key factors regulating reionization. Large simula- 
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Fig. 6. — Mass range of dark matter halos that can 
host galaxies. The curve M{z, f) gives the mass range 
M > M(z, /) that accounts for the fraction / of the 
expected total luminosity at redshift z. The five solid 
curves, from top to bottom, correspond to / = 0.2, 0.4, 
0.6, 0.8, and 1.0. If 30 to 50 particles are required to 
identify a halo in an N-body simulation, then the red 
and blue bands specify the particle mass resolutions that 
are needed to account for 50% and 100% the total lu- 
minosity, respectively. 

tion volumes are necessary to have a fair representation 
of the distribution of galaxies and ionized bubbles, and 
to find rare quasars. 

For studying the large-scale properties of reionization, 
the minimum simulation box size should be approxi- 
mately 100 Mpc on a side for two major reasons. 
First, as we discussed in the mean free path for 
ionizing photons is expected to be several tens of Mpc 
at 2 ^ 6 — 10. Thus, the simulation box size needs 
to be many times larger in order to have a fair sam- 
pling of the distributio n of ionized bubbles. Second, 
Barkana &: Loebl (|2004l ) showed that a box size of about 
100 Mpc is necessary to have a fair sample of dark 
matter halos where sources reside. Within this volume, 
the statistical fluctuations in the abundance of galaxies 
are small enough that is representative of the universe 
on average. Therefore, this volume is reionized no earlier 
or later than the universe on average by any significant 
amount. 

Quasars are much rarer than galaxies, especially at 
early times. The quasars observed at z ~ 6 in the SDSS 
have a very low number dens ity, roughly one f or every 
(1000 h^^ Mpc)3 of space fe.g. lFan et al.ll2006al ). Large- 
scale cosmological N-body simulations are able to find 



10 



massive dark matt er halos with m ass ~ 10"'^^ Mq, which 
host these quasars (jLi et al If the upper Hmit on 

the required box size is taken to be 1000 /i"^ Mpc, then 
the total mass within this volume, ^ 10^° Mq or about 
10^ times the mass of all matter in our galaxy, sets the 
corresponding mass limit. 

On small-scales, the situation is rather complicated. 
Fig. [S] shows the mass range of dark matter halos that 
is expected to host galaxies. For the redshift range of 
interest, halos with mass scale ~ 10^ M© h^^ and co- 
moving length scale ~ 10 h^^kpc are thought to be able 
to form stars more efficiently. Within these halos, the 
gas reaches an important temperature scale ~ 10^ K 
at which it can cool rapidly through atomic transitions. 
The dissipation of energy enables further gravitational 
collapse to much higher densities, allowing the forma- 
tion of giant molecular clouds and in these the forma- 
tion of stars. This sets the upper limit for what the 
minimum spatial and mass resolutions should be in or- 
der to locate where galaxies reside. The resolution has 
to be much better in order to resolve the gas structure 
within these halos. For an extreme example, we have to 
resolve down to tens of proper parsecs to directly form 
giant molecular clouds of order a thousand solar masses 
in gas. 

It is difficult to set resolution limits for radiation 
sinks because we do n ot fully understand wha t the 



absorbing systems are. iKohler fc GnedinI (|2007l ) have 
found that Lyman limit systems at z = 4 are gen- 
erally found near galaxies. They are typically a kpc 
in proper size and have densities roughly a thousand 
times larger than the universal average. While this 
gives us a sense of the required scales, the situation at 
higher redshifts could be different. For example, con- 
sider halos that were not massive enough to form stars 
efficiently. These mini-halos were very abundant com- 
pared to the source halos and their reservoir of absorb- 
ing gas could make t hem significant radi ation sinks (e.g. 
Haiman etallbOOll : [Shapiro et al.ll2004) . 

Having some idea of the physical scales, we consider 
the demand in dynamic range and computational re- 
sources. The extreme scenario, simulating the formation 
of giant molecular clouds within large cosmic volumes 
where rare quasars can be found, would require over 7 
orders of magnitude in length and 16 orders of magni- 
tude in mass. While adaptive simulations can achieve 
this dynamic range for any given galaxy, it is still far 
out of reach to accomplish this for every galaxy within 
a very large cosmic volume. 

For a more practical case, we consider simulations that 
capture the formation of the large-scale structure and 
the reionization of the intergalactic medium within a 
simulation box size of 100 h^^ Mpc. Given the distribu- 
tion of halos that host galaxies, radiation sources and 
sinks can then be modeled approximately by populat- 
ing these halos with assumed distributions for proper- 
ties such as source luminosities and sink opacities. For 



example, the source luminosity can be assumed to be 
proportional to the halo mass. In N-body simulations, 
approximately 30 to 50 particles are needed to properly 
identify a dark matter halo and measure its mass. To 
form halos with M = 10^ Mq then requires approx- 
imately 20 to 35 billion particles in total within the spec- 
ified volume. Only recently has this been achieved by 
high-resolution N-body simulations utilizing 512 to 2048 
processors, 2 to 4 trillion bytes (terabytes) of memory, 
and 100 to 200 tho usand c pu-hours on modern super- 
comp uters (|Shin et al..,2008i;llliev et ahlbOOStlTrac et al 
20081) . 



Ideally, an equal number of hydrodynamic elements, 
either Eulerian grid cells or SPH particles, would be used 
to represent the cosmic gas. Hydro simulations have not 
yet been run at this scale because they generally require 
at least twice as much computer memory and take over 
ten times longer to run compared to N-body simulations. 
The computational costs will also increase substantially 
with the addition of radiative transfer. This major mile- 
stone will be achievable within the next few years using 
supercomputers with of order 10,000 processors and tens 
of terabytes of memory. 

5. Computer Simulations 

The primary challenge for modern simulations is to 
address the open question of how the distribution and 
properties of sources and sinks affect the reionization 
of the universe. As previously discussed, it is not yet 
possible to probe the full range of relevant scales with 
any single simulation. High-resolution simulations with 
small boxes can resolve radiation sinks such as Lyman 
limit systems, but the number of sources is not large 
enough to be representative of the actual galaxy dis- 
tribution. On the contrary, large-box simulations have 
sufficient volume to be considered representative of the 
homogeneous and isotropic universe during the epoch of 
reionization, but generally have insufficient resolution to 
provide information on small-scale structure. 

For clarity of discussion, we classify reionization sim- 
ulations into two major categories: small-scale simula- 
tions that attempt to resolve known radiation sinks and 
large-scale simulations that attempt to account for the 
abundance of expected sources. In the first category, hy- 
dro + radiative transfer simulations directly model the 
evolution of the dark matter, cosmic gas, and radiation. 
These simulations resolve down to kpc scales or even 
smaller, allowing the modeling of two important phys- 
ical processes. First, high-density, cooHng gas can be 
found and turned into stars using simple, but physically 
plaus ible prescriptions (e.g. Katj 1992t ICen fc Ostriker 
19921 ). Second, the presence of high-density absorbing 
gas will act as Lyman limit systems. 

Simulations falling into the second category share the 
common traits of having large volumes and high source 
counts. Hydro -I- radiative transfer simulations are pre- 
dominantly not used, but instead hybrid approaches are 



11 



Fig. 7. — Visualization of the reionization process in small box simulations of iGnedin fc Fan (|2006l) . Opaque brown 
material represents neutral hydrogen. Yellow points are galaxies. The glowing blue color shows recently ionized gas. 
The three panels show three different moments: z = 8.1 (before the overlap of ionized bubbles), z = 6.3 (beginning 
of overlap), and z = 5.5 (after the overlap). 



taken. For example, high-resolution cosmological sim- 
ulations are first run to model density fields and halo 
distributions. Radiative transfer calculations are then 
performed by post-processing the density fields using 
sources modeled from the halos. Various hybrid tech- 
niques have been developed, including those in which 
unresolved small-scale physics can be incorporated in 
approximate ways. 

Studies of reionization necessarily overlap with many 
other directions of astrophysical and cosmological re- 
search. Constrained by the limited size of this review, 
we restricted our focus to computer simulations that 
model the whole reionization process and for which 
reionization is the prime subject of study. Thus, we do 
not discuss a large body of ground-breaking numerical 
work on modelin g the formation of the very first stars in 
the universe (see Bromm fc Larsonl 2004 ; NormanI 20081 
for recent reviews), and on simulating their effect on 
cosmic structure formation on small scales. Nor do we 
discuss an even larger body of numerical work on mod- 
eling the universe after reionization (see MeiksinlboOTl 
for a recent review). All of these simulations play an 
important role in developing our understanding of the 
high-redshift universe. 

5.1. Small-scale simulations that resolve sinks 

Small box simulations, that aim at resolving all 
sinks for ionizing photons, have been historically 
the first attempt to model the p rocess of reioniza- 
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incarnation of these simulations performed by one of us 
(NG) uses a box size of 8/i~^ comoving Mpc and reach 

pc in comoving 



units, or some 100 pc at 8 ^ z ^ 10. These simulations 
reproduce well the observed properties of L yman limit 
systems at z ~ 4 ( Kohler fc Gnedin 20071 ). However, 
one can never be sure that some other type of sinks 
exists during the reionization epoch; thus, the "small- 
box" simulations attempt to resolve all sinks, but there 
is no guarantee that they actually do that. 

An example of the small box simulations is shown in 
Fig. pfPl . The fact that the volume of these simulations 
is too small is immediately apparent from the figure: 
most of reionization is done by just two large ionized 
bubbles around two typical large galaxies (this is espe- 
cially apparent in the first panel) . The rest of galaxies in 
the simulation volume are small, dwarf galaxies; while 
they create a larger number of smaller ionized bubbles 
around them, these bubbles do not contribute signifi- 
cantly to the overall reionization process. 

Obviously, two is not a statistically significant sam- 
ple. In order to partially compensate for the small box 
of these simulations, the simulation volume is usually 
chosen not at random, but to be as close to an aver- 
age place in the universe as possible - such a choice 
requires creating several hundreds of initial conditions 
and choosing "the best one" for the actual simulation. 
Because of such special selection, these simulations can 
be used to reprodu ce some of the obse rvational data 
from SDSS quasars ( Gnedin fc "FM]l2006h that describe 
the average properties of the universe; for example, the 
average fraction of quasar light transmitted through the 
intergalactic medium, the mean fraction of neutral hy- 
drogen of the universe, etc. However, these simulations 
fail when they are used beyond computing simple aver- 
age properties. 

Another serious limitation of small box simulations is 



spatial resolution of better than 700/i ^ ^ 



^More still images and several animations are available at 
http: //home . f nal .gov/'gnedin/GALLERY/rei.p.html 
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Fig. 8. — A comparison of spectra of two z — 6.3 
quasars: a real observed spectrum of the SDSS quasar 
J1030+0524 and a synthetic spectrum from a large-scale 
simulation with "clumping factors" approach. Guessing 
which is which we leave to the reader as a practical ex- 
ercise. 

that reionization in them always proceeds in the "abun- 
dant sources" scenario, simply because the simulation 
volume is smaller than the mean free path due to Ly- 
man limit systems. Ionized bubbles, then, overlap be- 
fore they reach sizes comparable to the mean free path 
and never have a chance to stall, as in the "rare sources" 
scenario for reionization. 

Numerical study of reionization began with small box 
simulations, but by now they have largely completed 
their role. There is still one application remaining for 
this type of simulations, though, and that is their use 
in combination with very large volume simulations that 
utilize the "clumping factors" approach, which we dis- 
cuss next. 

5.1.1. Clumping factors approach 

Because the whole dynamic range, from Lyman limit 
systems to the scale of cosmic horizon, is currently 
unachievable in a direct numerical simulation, several 
approximate techniques have been developed recently. 
While we do not overview all approximate techniques 
here due to space limitations and our focus on simula- 
tions, one of these techniques - the so-called "clumping 
factors" approach - still falls into a general simulation 
category. In this approach the small-scales simulations 
that resolve Lyman limit systems are combined with 
large-scale simulations that include a representative vol- 
ume of the universe (we discuss such simulations in the 
next sub-section) , or even a larger volume up to several 
Gpc on a side to include the most rare quasars similar 
to ones observed by the SDSS quasar survey. 

The idea of clumping factors is simple. In a large-scale 
simulation of a Gpc-sized volume the spatial resolution 
is limited to a few Mpc at best, so a single resolution 
element (a mesh cell or a particle) in such a simulation 
is a uniform region of a few Mpc on a side. Certainly, 



the universe cannot be assumed to be uniform on this 
scale; a "clumping factor" approach uses the small-scale 
simulations to statistically describe the property of the 
universe on a few Mpc scale, and then this statistical 
description is used in a large-scale simulation as an ap- 
proximation to the correct mathematical term in the 
evolution equation. 

So far, the "clu mping factors" appro a ch has not been 
widely used yet (jKohler et al.l 120071 : iMcQuinn et al 



2007bh . Fig.[8]demonstrates the ability of this approach 
to reproduce the observed spectra of high redshift SDSS 
quasars, but only future work will demonstrate whether 
the ability of "clumping factors" aproach to model ex- 
tremely large dynamic range (more than 1,000,000) jus- 
tifies its approximate nature and a substantial compu- 
tational expense. 

5.2. Large-scale simulations that account for 
the abundance of sources 

Over this past decade, significant progress has been 
made towards understanding how the large-scale dis- 
tribution of galaxies affects the reionization of the 
universe. Early sim ulations (e.g. Ciardi et al. 20031 : 
Sokasian et al. I l2003h were restricted to small box sizes 
of 10 — 20 Mpc and small numbers of sources be- 
cause of the limited computing power then. In addition, 
the radiative transfer calculations were done in post- 
processing to simplify the computation. Despite the 
limitations, these simulations showed that the abun- 
dance and luminosities of galaxies clearly affected the 
reionization process and the development of ionized 
regions. The interesting results paved the way for con- 
tinuing studies. 

Since N-body simulations are less costly to run than 
hydro, they have been predominantly used in the last 
few years to generate density fields and halo distri- 
butions for hybrid modeling. The box size milestone 
of 100 Mpc f or reionization simulations has now 
been reached (e g. Illiev et al.l l2006bl IZahn et all 12007 : 
McQuinn et al. 2007a[ l. but onlv recently have both 



the box size and m ass resolution r e quirements been 
met simultaneou sly (|Shin et al.ll2008l: Illiev et all [20081: 
Trac et all l2008f) . However, the number of radiative 
transfer resolution elements used is considerably smaller 
because of the high computational cost. 

Fig. [9] is a samp le visualiza t ion fr om a modern large- 
scale simulation. Trac et al. ( 20081) used a hybrid ap- 
proach in which a high-resolution N-body simulation 
was first run to model sources and then hydro -I- ra- 
diative transfer simulations were run incorporating the 
sources. In the simulations, the radiative transfer of the 
ionizing photons proceeded such that large-scale, over- 
dense regions near sources are generally photo-ionized 
and photo-heated earlier than large-scale, underdense 
regions far from sources. The inhomogeneous process 
changed the thermal and ionization conditions, convert- 
ing the cold and neutral gas into a warm and highly 
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Fig. 9. — Visualization of the reionization process in a hydro + radiative transfer simulation from lTrac et alJ (j2008l ). 
The first four panels show the evolution of the ionized hydrogen density in a slice of size (100 Mpc)^ when the 
simulation volume is 25, 50, 75, and 100 per cent ionized. The fifth panel shows the temperature at the end of reion- 
ization while the last panel shows the redshift at which gas elements get reionized. Higher-density regions tracing the 
large-scale structure are generally reionized earlier than lower-density regions far from sources. At the end of reion- 
ization, regions more recently photo-ionized and photo-heated are typically hotter because they have not yet had time 
to cool. Higher resolution images and movies are available at ^http://www. cfa.harvard.edu/~htrac/Reionization. html 



ionized one, typically with temperatures ~ 10** K and 
neutral fractions ^ lO^"'. 

5.2.1. Ionized bubbles 

The large-scale distribution and properties of ionized 
bubbles have been studied in more detail in the last 
few years. Recent simulations (e.g. Iliev et ahl 2006b ; 
Zahn et al.ll20 07: Mc Quinn et al. l2007bHLee et al.ll200g ; 
Shin et al.l f200& ,Croft fc Altavl boOSi ) support the ba- 



sic picture for the development of ionized bubbles de- 
scribed in US] Ionized bubbles originate within halos 
hosting sources and expand outwards, such that higher- 
density regions near sources are ionized earlier than 
lower-density regions far from sources. Bubbles will 
merge with one another until they overlap and fill all 
of space. 



McQuinn etahl (|2007bf) and ICroft fc Altavl (|2008l ) 

have simulated a range of models for sources and sinks 
and are in agreement that sources more strongly influ- 
ence the development of ionized bubbles. The properties 
of sources can affect the morphology in several major 
ways. Ionized bubbles are generally aspherical, as seen 
in Fig. 121 but rarer sources tend to generate larger and 
more spherical bubbles. Being more strongly clustered, 
rare sources have many other smaller sources nearby 
that contribute photons and help ionize a larger vol- 
ume. On the other hand, lower luminosity sources can 
significantly add to the abundance of small bubbles, 
if they are not already embedded within larger bub- 
bles. The lower limit on luminosity will depend on the 
efficiency of star formation in smaller mass halos. 

While not the dominant factor, the distribution of 
sinks can still have strong effects on reionization. The 
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Fig. 10. — Size distribution of ionized bubbles when 
reionization is half completed. The results, from two 
ray-tracing codes (McQuinn et al. 2007b in red; Trac & 
Cen 2007 in green) and a semi-analytic code (Zahn et 
al. 2007 in blue), are in good agreement overall. They 
have similar characteristic peak sizes, but some differ- 
ences are seen at small and large radii. Figure is cour- 
tesy of O. Zahn. 

abundance of sinks, from Lyman limit systems and mini- 
halos, restricts the mean free path of ionizing photons. 
Once the mean free path is shorter than the bubble 
size, the ionized regions stop growing, thus delaying the 
reionization process. In the presence of stronger sinks, 
the distribution of bubble sizes shifts towards smaller 
bubbles. 

For a partially ionized universe, the probability dis- 
tribution or histogram of bubble size generally increases 
towards smaller sizes. There are relatively more small 
bubbles than large bubbles, reflective of there being 
more low-mass halos than high-mass halos. This is 
the nature of the hierarchical universe in which larger 
and rarer objects are assembled from smaller and more 
abundant objects. Note that there is not a one-to-one 
correspondence between bubbles and sources since ion- 
ized regions can have one or more galaxies embedded 
within. The bubble size distribution can also be cal- 
culated by weighting each bubble by its volume rather 
than a simple count. In this version, the distribution 
has a characteristic peak size that shifts toward larger 
scales as reionization progr esses and ionized bubbles 



('F urlanetto et al.H2004h 



Fig. [TUl shows a sample result on the bubble size distri- 
bution from a code coinparison project. Two ray- tracing 
codes (|McQuinn et al.ll2007bl : iTrac fc Cenll2007l) and a 



semi-analytic code (jZahn et al.l 120071 ) have calculated 
the reionization process, starting from identical initial 
conditions of a 100 h"^ Mpc box. The results shown are 
taken at the same redshift, when each simulation has 
reionized half the volume. The three calculations arc in 
good agreement overall, including having the same char- 
acteristic peak size, but small differences are present due 
to the different treatments for radiative transfer. 

In the last few years, semi-analytic models for reion- 
ization (e.g. Barkana fc Loebl 2004t Furlanetto et al 



2004 IZahn et al.. .2007.: .Mesinger fc Furlanettol 120071: 



Choudhurv et al. 20091 ) are predominantly b ased on a 



technique called the 'excursion set formalism' l|Bond et al 
In these schemes, the source distributions and 
ionization fields are approximately derived from the 
density fields. Basically, sources are first located within 
the density field by asking whether there is enough 
surrounding mass at a given point in space. Ionized 
bubbles are then found by asking whether a given re- 
gion has enough sources within it to ionize the con- 
tained mass. Since the propagation of photons is not 
directly followed, the algorithms tend to be much faster 
and use less memory than radiative transfer algorithms. 
Semi-analytic models have been demonstrated to be in 
good agreement with radiative transfer simulations, but 
there still are differences to be worked out. They are 
very useful for studying reionization, especially when a 
large number of models needs to be run. 

5.2.2. 21 cm radiation 

Neutral hydrogen from the epoch of reionization can 
be detected through high-sensitivity radio observations. 
In a hydrogen atom, the proton and electron both have 
spins and the alignment of the spins can flip between be- 
ing parallel and antiparallel. This spin flip, sometimes 
referred to as a hyperfine transition, has an energy level 
difference corresponding to a rest wavelength of 21 cm. 
During reionization, this signal is seen in emission or 
absorption relative to the CMB and the brightness of 
the signal can be used to probe the characteristics of 
neutral regions and the delineation with ionized regions. 
Studying the epoch of reionization and the high redshift 
universe with 21cm radiation i s a frontier topic i n cos - 
mol( 3gy. We refe r the reader to Furlanetto et al. ( 20061 ) 
and iLoebl |2008h for recei it reviews. 

Recent simulatio ns fe.g.'Mellema et al.^ 2006bt Lidz et al 

20071: ISantos et al l boOS : Back et al 2009|) have studied 

the contrast between neutral and ionized regions by 
calculating the expected 21cm brightness, usually ex- 
pressed as a "brightness temperature" , the quantity 
widely used in Radio Astronomy. The brightness tem- 
perature depends on the neutral hydrogen density and 
the gas spin temperature, which characterizes the frac- 
tion of parallel to antiparallel spin states. The signal is 
seen in emission when the spin temperature exceeds the 
CMB temperature and in absorption if the opposite is 
true. 
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Fig. 11. — Maps of the 21 cm brightness temperature (mK) showing the transition frq r n abso rption (blue) to emission 
(red) in a shoe of size (100 ft-^^Mpc)^, taken from the simulation of Santos et al. ( 20081) . The first panel shows 
the absorption signal from a very early stage of reionization when the spin temperature was less than the CMB 
temperature. In the second panel, neutral regions, surrounding ionized bubbles, that have been heated by high- 
energy X-ray photons have switched from absorption to emission. The last panel shows the emission signal when 
the simulation volume is half ionized, but the spin temperature exceeds the CMB temperature everywhere. Figure 
is courtesy of Alexandre Amblard and Mario Santos. 



The evolution of the gas spin temperature depends on 
three mechanisms: collisions within the gas, scattering 
with the CMB, and scattering with Lyman alpha pho- 
tons. The latter depend ence is commonly known as th e 
Wouthysen-Field effect (' Wouthuvse^ll952l: lFieldlll959l) . 
More recent simulations have made the effort to calcu- 



late th e evolution of the spin temperature. ISantos et al.l 
lioOS*) calculated the rise of the spin temperature by 
including heating of neutral gas by high-energy X-rays 
and pumping of s pin states by Lyma n alpha photons, as 
sources turn on. Baek et al. I (|2009l) performed a more 
exact calculation of the radiative transfer of Lyman al- 
pha photons, taking into account the resonant scattering 
with neutral hydrogen. These inclusions are important 
in earlier stages of reionization when sources are few. 

Fig, [m shows the transitioning of the 21 cm brightness 
temper ature from absorpti on to emission in the simula- 



tion of ISantos et al 



( 2008f ). Initially, the signal is seen 
in absorption because the neutral hydrogen gas is cold 
and the spin temperature is less than the CMB tem- 
perature. As sources turn on and produce radiation, 
the signal vanishes in the ionized regions. The X-rays 
and Lyman alpha photons will travel beyond the ionized 
regions and raise the spin temperature of surrounding 
neutral regions, resulting in the transition from absorp- 
tion to emission. As reionization progesses, the emission 
regions will shrink in size. 

The primary statistic for studying the 21 cm signal 
has been the power spectrum of the brightness temper- 
ature field. The power spectrum quantifies how fluc- 
tuations in the brightness field are correlated with one 
another. The power spectrum of a field 5(x) is normally 
obtained by first Fourier transforming the field and then 
calculating the average power (|i5(k)p) for modes with 



wavenumber fc = |k|. It is also the Fourier transform of 
the two-point correlation function, which was discussed 
in S}3l and thus, a measure of spatial correlations. While 
the 21 cm signal comes from the neutral regions, the rel- 
ative fluctuations in the brightness field are commonly 
discussed in the context of ionized bubbles. 

The 21cm power spectrum has some interesting fea- 
tures. On scales large compared to the characteris- 
tic bubble size, the 21cm power spectrum is approxi- 
mately proportional to the power spectrum of the mat- 
ter. This correlation is due to the fact that ionized 
bubbles approximately trace the large-scale structure, 
as seen in Fig. [9l The scaling or bias factor comes 
from the fact that bubbles can cluster differently than 
the matter. The power spectrum also has a charac- 
teristic peak scale corresponding to the characteristic 
bubble size. As reionization progresses and the bub- 
bles merge, the peak shifts toward larger scales and the 
power spectrum changes shape. Several groups of ra- 
dio astronomers are actively developing instruments to 
measure this signal early in the next decade. These ob- 
servations will provide a wealth of information on the 



great detail (e.g Furlanetto et al.|l2006 


; Lidz et al. 20081 


Pritchard & Loeb 20081 lBarkanall2008 


)■ 



6. Conclusions 

Computer simulations of reionization have progres- 
sively achieved better physical realism and dynamic 
range over the past decade. The constant improvement 
of numerical algorithms combined with the rapid ad- 
vancement in supercomputing resources have fueled the 
growth. With N-body and hydro simulations, we can 
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now robustly evolve the dark matter and cosmic gas to 
form the large-scale structure of the universe. In con- 
junction, radiative transfer algorithms, based on mo- 
ments, Monte Carlo, and ray-tracing methods, have en- 
abled us to accurately solve the propagation of ionizing 
photons. While various implementations exist, good al- 
gorithms possess the important features of conserving 
photons and scaling linearly with the number of resolu- 
tion elements. 

Armed with modern tools, we have made significant 
strides in understanding how the universe was reionized 
and how sources and sinks influenced the process. The 
development of ionized bubbles, first around sources, 
followed by expansion to merge with other bubbles, and 
finally overlapping all of space, is now a quantifiable 
theory. Statistical tools have been developed to charac- 
terize the morphology and clustering of bubbles, and to 
discriminate between the possible scenarios for reioniza- 
tion. 

One of the major challenges remaining is to develop 
better physical models for sources and sinks that can 
be tested against high-redshift observations. Ideally, 
our ignorance of complex astrophysical processes can 
be represented by a small number of unknown, but con- 
strainable, parameters. The task of a computational 
cosmologist then becomes to ensure the high fidelity of 
the numerical solution given these parameters, so that 
modeling the process of reionization reduces to search- 
ing the parameter space. Since reionization is a complex 
process with no general analytical solution, the check 
on simulations requires having numerical codes converge 
with one another. 

Wo will have better observational constraints on reion- 
ization from ongoing surveys of high-redshift quasars 
and young, star-forming galaxies. The latter is espe- 
cially important for learning about sources. Upcoming 
radio observations of the 21 cm radiation from neutral 
hydrogen have the potential to be the best probe of the 
epoch of reionization. In addition, higher-resolution ex- 
periments at microwave wavelengths will be able to mea- 
sure the fluctuations in the CMB temperature caused 
by scattering with electrons produced by the reioniza- 
tion process. A wealth of information is expected over 
the next decade. Continuing progress on the theoretical 
front is therefore necessary to compliment the exciting 
promise of future observations. 

We thank Alexandre Amblard, Xiaohui Fan, Ilian 
Iliev, Mario Santos, Volker Springel, and Oliver Zahn 
for kind permission to use their figures in this review. 
H.T. is supported by an Institute for Theory and Com- 
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